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Abstract 

We show how to obtain a fast component-by-component construction algorithm for higher 
order polynomial lattice rules. Such rules are useful for multivariate quadrature of high- 
dimensional smooth functions over the unit cube as they achieve the near optimal order 
of convergence. The main problem addressed in this paper is to find an efficient way of 
computing the worst-case error. A general algorithm is presented and explicit expressions 
for base 2 are given. To obtain an efficient component-by-component construction algorithm 
we exploit the structure of the underlying cyclic group. 

We compare our new higher order multivariate quadrature rules to existing quadrature 
rules based on higher order digital nets by computing their worst-case error. These numer- 
ical results show that the higher order polynomial lattice rules improve upon the known 
constructions of quasi-Monte Carlo rules based on higher order digital nets. 
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1 Introduction 



In this paper we are concerned with quasi-Monte Carlo rules, which are equal weight multivariate 
quadrature rules (or cubature rules) 

N-l 

used to approximate multivariate integrals over the s-dimensional unit cube 

/(/):= / f{x)dx. (2) 

^[0,1] = 

In contrast to the Monte Carlo method, which samples the function / randomly in its domain, 
the integration nodes {xh}^ZQ used in the quasi-Monte Carlo rule Q are chosen deterministically. 
The convergence of the error of the Monte Carlo method is 0{N~^^'^) (in distribution), while 
the worst-case error for quasi-Monte Carlo is 0{N~°'{\ogNY'^) [4, 5, 11, 15] and the random- 
case error for randomized quasi- Monte Carlo is 0(iV-"-i/2(iog Ar)^("+i)) [3^ 22]. The latter two 
methods require that the integrand has smoothness a > 1 (which means for instance that the 
integrand has square intcgrablc partial mixed derivatives up to order a in each variable) , whereas 
the Monte Carlo method requires only that the integrands have finite variance. 

Different types of point sets for quasi-Monte Carlo rules exist. Those of interest here are digital 
nets [11, 15]. These can be divided into (classical) digital nets [15], which achieve a convergence 
rate of 0{N~^ {log NY) [15] for integrands of bounded variation, and their extension called higher 
order digital nets [4, 5], which achieve a convergence rate of 0(A''~"(log A^)"*) for integrands 
which have square integrable partial mixed derivatives of order a > 1. In [5, Section 4.4] an 
explicit method for constructing such higher order digital nets, based on a classical digital net, 
can be found. The method in the current paper gives an alternative construction for a higher 
order version of a specific type of digital net, namely polynomial lattice point sets, see [15, 
Section 4.4] or [11. Chapter 10]. Constructions of classical polynomial lattice point sets based on 
a worst-case error criterion have previously been studied in [8] . Note that we call a quasi-Monte 
Carlo rule whose underlying quadrature points are polynomial lattice points a polynomial lattice 
rule. 

Polynomial lattice point sets were generalized in [10] to obtain higher order polynomial lattice 
point sets. In [7] existence results on higher order polynomial lattice point sets were compared 
to the explicit construction of higher order digital nets [5] in terms of their t-value (a certain 
quality measure). For some values of dimension s and/or smoothness parameter a the higher 
order polynomial lattice point sets have a better existence bound than the best results which 
can currently be obtained using the explicit construction of higher order digital nets from [6]. 
The same is also true for classical digital nets, sec [14, 24]. These findings motivated the quest 
for an explicit construction of higher order polynomial lattice rules in [2]. The construction 
employed there is an algorithm originally proposed for the construction of (integer) lattice rules, 
namely the component-by- component construction algorithm, see, e.g., [12. 13, 25]. The higher 
order polynomial lattice rules so constructed achieve nearly optimal rates of convergence. For 
analogous results on polynomial lattice point sets see [8] (and also [9] for more background). 

Straightforward implementation of the component-by-component (CBC) algorithm is however 
very costly with respect to computational time, hence methods for reducing the computational 
cost are needed. The fast component- by-component algorithm, introduced in [20], uses fast 
Fourier transforms (FFTs) to speed up the calculations. Some notes concerning the application 
of the fast algorithm to the construction of polynomial lattice rules were already made in [21], 
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with a more detailed analysis in [19]; see also [11, Section 10.3]. In this paper we will adapt the 
fast algorithm for higher order polynomial lattice rules. To do so, we find a closed form for the 
worst-case error of our function space (where we consider the worst-case error as a function of 
the quadrature points). We show that our algorithm has a computational cost of 0{sN°'a\ogN) 
using 0{N°') memory, compared to 0(s^iV""'"^) for the straightforward implementation of the 
algorithm in [2], using the same amount of memory. This speedup makes it possible to obtain 
higher order polynomial lattice rules for moderate dimensions and numbers of points. In the 
section on numerical results we provide constructions of higher order polynomial lattice rules in 
base 6 = 2 up to dimension 10 and up to 4096 points. These numbers could be increased with 
more computational effort, but we have to remark that the search space grows exponentially 
with respect to the smoothness parameter a. 

The efficient calculation of the worst-case error of our function space is an essential ingredient 
in such an algorithm. We show that the kernel function associated with the worst-case error can 
be evaluated at a point x in time 0{an), where a is the smoothness of the space and a; is a 
rational number < v < 6". Moreover, in the case of the greatest practical importance, 

i.e., where the base equals 2, we show explicit expressions for smoothness 2 and 3 which are exact 
for any real x G [0, 1) (see Corollary 1). 

We compare the performance of higher order polynomial lattice rules constructed using our 
fast component-by-component algorithm to the explicit construction as outlined in [5] and find 
that the new algorithm performs better in the cases considered. Finally, for the benefit of the 
reader, we present some limited tables of higher order polynomial lattice rules constructed using 
the fast component-by-component algorithm, allowing the reader to apply the rules presented in 
this paper to problems of interest and to verify implementations of the algorithm. 

In the next section we provide the reader with some background, and notation, on Walsh 
spaces, digital nets, and the worst-case error. More detailed information can be found in [11] 
and [2] , where also bounds on the worst-case error for higher order polynomial lattice rules were 
proven. In Section 3 we show how to efficiently calculate the worst-case error and how the 
construction of higher order polynomial lattice rules can be done using the fast component-by- 
component approach of [19, 20]. 

2 Background 

We first introduce some notation. Let No = {0, 1,2,.. .} denote the set of non-negative integers 
and N = {1,2,3,...} the set of positive integers. Further we need to be able to consider a 
non-negative integer A: S No in its unique base b representation: 



where G {0, . . . , 6 — 1} are the base b digits of k and 7^ 0; a = for a; = 0. Note that the 
base, 6, is considered a fixed integer throughout. Moreover, in the further development in this 
paper, b will be prime. We will be specifically interested in the non-zero base b digits of k. The 
number of non-zero base b digits of an integer k will be denoted by ^fc; where ^0 = 0. We can 
then represent fc G No uniquely as 



a 




(3) 



#fe 




(4) 
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where now G {1,...,6 — 1} and we demand ai > • ■ • > a^k ^ 0. Thus is the most 
significant base h digit of k. For real x G [0, 1) we write its base 6 representation 

CO 

X={0.U2■.■)b^Y.^^h-\ (5) 

i=l 

where e {0, . . . , 6 — 1}. This representation is unique in the sense that we do not allow an 
infinite repetition of the digit — 1 to the right. 



2.1 A function space based on Walsh series 

For fc G No the one-dimensional A;th Walsh function in base 6, walfc : [0, 1) — !■ C, is defined by 

walA;(a;) := exp(27ri (^iKq H h ia+iHa)/b), (6) 

where we have used the base h digits of x and k as given in (3) and (5). Note that Walsh functions 
(in base h) are piecewise constant functions. For dimensions s > 2 and vectors k — (fci, . . . , fcg) G 
Ng and a; = (cci, . . . , Xs) G [0, 1)'' we define walfc : [0, 1)'' -> C as 

s 

walfe(a;) := ]^ walfc. (xj). 

The integrand functions in this paper are assumed to have an absolutely convergent Walsh 
series representation 



where the Walsh coefficients are given by 



fk^ /(a;) walfc (a; ) da;. 



'[0,1 

Note that the Walsh functions form a complete orthonormal system of L2{[0,lY)- For more 
information on Walsh functions and their properties we refer to [11, Chapter 14 and Appendix A]. 

In the following wc define a function space by demanding a certain decay rate of the Walsh 
coefficients. To do so, we introduce some further notation. We define, for a fixed integer a > 1 
and a fixed sequence of positive weights 7 = {71, 727 • • •} (in the sense of [26]), 

rah,k):^l[j,r^ik,), r„(fc) (7) 

where for k = {ki, . . . , kg) we set u = {\ < j < s : kj ^ Qi\ and where we used the first a 
positions oi + 1, . . . , a^u + 1 of the non-zero base 6 digits of k with the notation defined in (4) 
in the one-dimensional definition of r^. For fc = = (0, . . . , 0) we set 

r„(0) = 1. 

We are now ready to specify which functions are in our function space VVq^s,^, which was also 
used in [2, 5]. For functions / G yVa.s,-^ wc define the norm 

[|/[[w„,.,. sup (8) 
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Then yVa,s,~f consists of all functions / € L2([0,l]^) for which ||/||wc,»,-, < oo. The Walsh 
coefficients of / G yVa,s,-y therefore satisfy a certain decay criterion, namely 

l7fc| < ||/||w„..,, r„(7,fe) VfceNg. (9) 

It is clear that larger values of a might increase the norm of a function /, i.e., ||./||wa 3 l£ 
||/||w^, ^ for a < a' . The weights 71, 72, . . . arc used to describe how anisotropic the space is. 
Usually it is assumed that 71 > 72 > ■ ■ • > 0, meaning that the first dimension is more important 
than the second one and so on. Under certain conditions on these weights, it can be shown that 
numerical integration is tractable in the number of dimensions, sec, e.g., [11, 16, 17]. 

It is of course important to have an understanding of which functions exactly arc in such 
a Walsh space Wa,s,-y with smoothness parameter a. This analysis has been done in [4, 5, 6]. 
Classically, one is interested in (smooth) functions / : [0, 1]** — >■ R for which all mixed partial 
derivatives up to order a in each variable are square integrable. This is a Sobolev space of 
smoothness a which is often considered for this type of problems. In [5, 6] a continuous embedding 
of certain Sobolev spaces into Wa,i5,-y was shown. Consequently, the results we are going to 
establish in the following for functions in Wa.s,-f also apply automatically to what we normally 
consider as "smooth" functions, for instance, functions which have square integrable partial mixed 
derivatives up to order a in each variable. One of the simplest type of functions in this space are 
multivariate polynomials which make up nice testing examples for computer implementations. 

2.2 Higher order digital nets 

Higher order digital nets were introduced in [5]. Higher order polynomial lattice point sets, which 
are the focal point of this paper, are a special class of higher order digital nets. For that reason 
and since we will compare the explicit construction for higher order digital nets from [5] with the 
construction given in this paper, we will review the necessary details here. For more information 
we refer to [11, Chapter 15]. 

For a prime number b we always identify Ft,, the finite field with b elements, with Zj, = 
{0, . . . , 6 — 1} endowed with the usual arithmetic operations modulo b. 

First we define higher order digital nets using the digital construction scheme. As we will need 
to be able to identify integers with vectors over a finite field by using its base b representation, 
and then later have to be able to consider vectors of integers as well, we will denote a vector over 
a finite field Ff, by h, in contrast to vectors over Z or R, which will be denoted by h. 

Definition 1 (Digital construction scheme of a digital net over Ff,). Let 6 be a prime and let 
n,m,s > 1 be integers, where n > m. Let Ci, . . . , Cs be n x m matrices over the finite field F;, 
of order b. Now we construct fe™ points in [0, 1)''': for < /i < 6™, identify each h = X^I^lo^ 
with a vector over the finite field 

h:= {ho,...,h^^,f e¥l^. 
For I < j < s multiply the matrix Cj by h using arithmetic over F^ to obtain a vector i/hj G FJ,': 

Cj h vhj = (y/ijM, . . . , Vhj.nV e f;', (10) 

from which the hth point of the digital net is found by interpreting the coordinates of ijhj as 
the base b digits of x^j'. 

n 
i=l 
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Now set Xh ~ {xh,i, ■ ■ ■ ,Xh,s)^ S [0, 1)^ to be the hth point. The set {xq, . . . ,a;fc™_i} is called 
a digital net over F;, with generating matrices Ci , . . . , Cs . 

This definition of a digital net generalizes the classical construction scheme, e.g., [15], on 
which classical digital (t, m, s)-nets are based upon, by allowing for generating matrices which 
are not necessarily square. The generating matrices Cj are of size n x m, and so, the number of 
rows n determines the resolution at which the points of the net are placed in the unit cube, i.e., 
all base b digits after position n are zero. The integration error then behaves like 0{b~'^k°"'), see 
[5], where a is the smoothness of the integrand and k is the strength of the net (in accordance 
with the respective property of classical nets). The strength of the net is defined via linear 
independence properties of the rows of the generating matrices, see [4, 5]. For higher order nets 
one can achieve k ~ min(am, n) and hence, provided that n > am, one obtains a convergence 
order of 6""™ (am)"* x A^^"(logiV)"'*, where = is the number of quadrature points. 

We now explain the explicit construction of a higher order digital net in s dimensions for 
a maximum smoothness d as described in [5]. The explicit construction starts from a given 
{t' , m, sd)-net in base b, that is, a classical digital net in sd dimensions for which the generating 
matrices Ci, . . . , Csd g F"^™ are known. 

From these sd given matrices, s new generating matrices cj'*'' are constructed of size dm x m 
by vertically stacking the first rows from the group of d consecutive matrices Cq-ij^j+i, . . . , Cjd, 
then the second rows of the same d matrices and so on, until all dm rows have been stacked. 
More precisely, let Cj ~ {cji . . . cj^)^, where cjf, denotes the fcth row of the matrix Cj. Then 



C] = (C(^_i)d+i^i . . . cj^,i . . . cj^_i)d+i,m ■ ■ ■ ^Jd,J^- For more information on these higher 
order digital nets we refer the reader to [5] and [11, Chapter 15]. 



An important concept for the error analysis in the next section is the dual net. It defines the 
set of Walsh coefficients which are not integrated exactly by the digital net and can therefore be 
used to write down the integration error. 

Definition 2 (Dual net). For a digital net over F(, with generating matrices Ci, . . . , Cs G FJ^^™ 
we define its dual net by 



where for a scalar component k = X^i^o '^i ™ ^ define an associated vector over the finite 



2.3 The worst-case error 

We define the worst-case error of numerical integration using a cubature rule Q for functions in 
a Banach space J- by 



We now assume the cubature rule Q to be a quasi- Monte Carlo rule (1) using a (higher order) 
digital net as its node set and denote it by Q""*. For any / having an absolutely convergent 
Walsh series representation we can write the integration error for Q"°' as a sum over the dual 
net to obtain 



I?(Ci,...,a) { 




field ^= (ko,...,k„)^ GF^. 



e(g,^):= sup |/(/)-Q(/) 



ll/lt^<i 




(11) 
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For / G Wa,s.-y wc can now use (9) to obtain 

|/(/)-0"^*(/)l<ll/l|w„_ J2 ^"('^'^)- (12) 

o=ikev 

Since we can obtain equality for a worst-case function ( G Wa.s,-y having Walsh series represen- 
tation 

C(^) = ''"(T'fc) walfc(a;), 

feGN= 

we find the following expression for the worst-case error in Wa.s,-y for a quasi-Monte Carlo rule 
(5"°' based on a higher order digital net in base b: 

e(Q"^*,>Vc«..,-y) = ^-h,k). (13) 

The cubature rules in this paper will be constructed in such a way that they have a worst-case 
error which is near optimal for the given function space VVa,s,-y. For a given value of a > 1 
the worst-case error behaves like 0(A^~"(logiV)"'*) for N integration nodes (see [5]) which is 
essentially best possible according to a lower bound from Sarygin [28]. 



2.4 Higher order polynomial lattice rules 

In [10] the classical polynomial lattice rules [15] were generalized to form higher order polynomial 
lattice rules. Just like classical polynomial lattice point sets are a special class of digital nets, 
higher order polynomial lattice point sets are a special class of higher order digital nets. For 
simplicity, we define the (higher order) polynomial lattice rules over a finite field Ff, of prime 
order b only. The main object for the construction of polynomial lattice rules are formal Laurent 
series, i.e., expressions of the form Wi X~'\ where £ G Z and m; G Ff,. Wc denote the set of 

formal Laurent scries by F(,((X~^)). These Laurent series then need to be mapped to integration 
nodes over the unit interval [0, 1). Define the map w„ : Vb{{X~^)) [0, 1) by 

n 

2— max(£, 1) 

Similar to the case for digital nets, we now need to be able to identify an integer h with a 
polynomial in Fb[X] by considering h in its base b representation, the associated polynomial will 
be denoted by h{X). The details are given in the following definition. 

Definition 3 (Polynomial lattice rule). Let b be prime and 1 < m < n. For a given dimension 
s > 1, choose p{X) G ¥b[X] with deg(p) = n > 1 and let gi(X), . . . , qs{X) G ¥b[X]. Now we 
construct 6™ points in [0, 1)": for < /i < 6™, identify each h = X^I^o^ ^i^h a polynomial 
over Ff, 

m — 1 

hiX) := h,X'e¥h[X]. 

i=0 

Then the /ith point is obtained by setting 

/ f hiX)q^iX) \ f hiX)qAX) \\ 

A quasi-Monte Carlo rule using this point set is called a polynomial lattice rule. 
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One obtains classical polynomial lattice rules from Definition 3 by taking n = m. For sim- 
plicity we will assume that p{X) is irreducible over Ff,, though this assumption could be removed 
by a more intricate analysis. We define 

Gb.n = MX) e ¥b[X] \ {0} : deg(z;) < n}, 

which will be the set from which we will select the generating polynomials qj(X). Clearly, as 
p{X) is irreducible and deg(p) = n, this equals the multiplicative group 

Gb,n = (F,[X]M^))x = {g{Xf : < /3 < 6" - f}, 

where g{X) is a generator for the multiplicative group {¥i,[X]/p{X))^ , e.g., we can take g{X) = 
X when p{X) is primitive. 

Since a polynomial lattice point set is a special case of a digital net, we can find the generating 
matrices Ci, . . . , G FJ^^™ from the generating vector q{X) = {qi{X), . . . , qs{X)). For 1 < j < 
s consider the Laurent expansions 

Iji^) _ „,(J) Y- 



PiX) 



f?) 

Then the elements cj:^ of the n x m generating matrix Cj over Ft, are given by 

for f < fc < n, < € < TO — 1; see, e.g., [11, Section 10.1]. 

In (13) we used the dual net to obtain the worst-case error. In the case of a polynomial lattice 
rule the dual is given in the next definition. (We use the convention deg(O) = — oo.) 

Definition 4 (Dual polynomial lattice). A polynomial lattice with generating vector q{X) = 
{qi{X), . . . , qs{X)) 6 (F;,[X])''' modulo p{X) e Wb[X] has a dual polynomial lattice 

V{q{X),p{X)) ;= J fe e Ng : ^ kj{X) qj{X) = a{X) (mod p{X)) with deg(a) < n - i 

A proof for the equivalence of Definition 2 and Definition 4 for polynomial lattices follows 
from [11, Lemma 15.25]. 

Specifically for a polynomial lattice rule with generating vector q{X) = (qi{X), . . . ,qs{X)) 
modulo p{X) having 6™ points, it follows from (13) that its worst-case error in Wa^s.-y satisfies 

efc™,„(q(X),p(X)) = ^ r„(7,fe), (16) 
with V the dual polynomial lattice. 



2.5 The component-by-component construction of higher order poly- 
nomial lattice rules 

The component-by-component construction algorithm was introduced by Korobov [12], see also 
[13, Theorem 18, p. 120], and later re-invented in [25] to construct the generating vector of an 
integer lattice rule. This algorithm first finds the optimal one-dimensional generating vector. 
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which is subsequently extended in an optimal way to a two-dimensional generating vector and so 
on. Algorithm 1 spells out the details for the construction in the case of higher order polynomial 
lattice rules. 



Algorithm 1 General form of CBC construction of higher order polynomial lattice rules 

Input: base b a prime, number of dimensions s, number of points 6™, smoothness a > 1, and 
weights 7 = (7j)j>i 

Output: Generating vector q{X) = {qi{X), . . . , qs{X)) £ Gl^ 

Choose an irreducible polynomial p{X) G Fh[Ar]. with deg(p) = n and n = am 
for d = Ito s do 

Set qdiX) £ Gb,n by minimizing a.{{qi{X), . . . , qd{X)),p(X)) as a function of qd{X) 
end for 

return q(X) = . . . , q,(X)) 



The analysis of the component-by-componcnt algorithm adjusted to the case of higher order 
polynomial lattice rule was done in [2]. The following theorem shows that Algorithm 1 achieves 
almost optimal rates of convergence. For a proof we refer to [2]. 

Theorem 1. Let b be prime, s,n G N and p{X) £ ¥b[X] be irreducible with deg(p) = n, a > \. 
Suppose {qi (X), . . . , qa{X)) £ ^ is constructed using Algorithm 1. Then for all d ~ 1, . . . , s 
we have a bound on the worst-case error as follows: 

e,™,„((qi(X), . . . , qd{X)),p{X)) < ^— n (l + 37j/'a,a,r)' , V 1 < r < a, 

where 



Q-i I a— 1 if T = I, 

■= n + { {b- mb - ir-^ - - 1)"-^) 

[ - l)"-i ^ ■ 

Formula (16) for the worst-case error is not in a usable form for computation due to the 
infinite sum. The next lemma shows how to obtain a closed-form expression which resembles 
the formula for the worst-case error as it appears when the space of integrands is a reproducing 
kernel Hilbert space, see [1]. 

Lemma 1. The worst-case integration error in Wa.s^-y^ ct > 1, associated with a polynomial 
lattice rule with generating vector q{X) ~ {qi(X), . . . ,qs{X)) modulo p{X) having 6™ points 
satisfies 



e,™.„((<7i(X),...,g,,(X)),p(^)) = -l + ^ ^ n(l + 7,w„(x,.,,)), (17) 

h=0 j=l 

where, using (7), 

oo 

,{x) := ^ra{k) walfc(a;). (18) 



k=l 
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Proof. Wc make use of the character property of digital nets (see [11, Lemma 4.75]). When 



{^h}^h=o ^ points of a digital net (or a polynomial lattice rule) and V is its dual net, then 

1 Jl iffcGP, 

— > walfc(a;/i) = < 
0™ f—' otherwise. 

h=0 K 

Thus, starting from (16), we obtain 

eb'^^a{{qi{X),...,qs{X)),p{X)) = ^ r„(7, A;)— ^ walfc(a?,,) 

= "1 + XI H ^"(7,fc) walfc(a;/,) 

h=0 fcGNg 

6"' - 1 s 

□ 

As in [20] we can now do a basic operation count for the computational cost of Algorithm 1. 
In comparison to [20], the analysis is a little bit more involved here as the evaluation of (17) 
involves calculating 

(^Vn ( ^^^(^)^^ ) ) for /i = 0, . . . , 6" - 1 and ah q{X) e Gf,,„, (19) 

in each iteration. Assuming c^^ = Cuj{a,n) to be the cost of evaluating uJa{x) and c„ = Cy(n) 
the cost of mapping and calculating the Laurent expansion as well as calculation the polynomial 
product, the cost of a straightforward implementation of Algorithm 1 is 0(s^&"&™(ctj + Ct,)) where 
n = am. However, calculating uJa{vn{h{X) q{X) /p{X))) efficiently for given h{X) and q{X) is 
an important issue which is solved in the next section and in practice it would be inefhcient to 
calculate these values on the fly whenever needed. To that end we model the algorithm into a 
more tangible form using the techniques from [20, 21] to obtain a fast componcnt-by-component 
algorithm that makes use of a circular convolution which can be calculated by means of fast 
Fourier transforms (FFTs). 



3 Fast construction of higher order polynomial lattice rules 

The exposition here mainly follows the techniques from [20, 21], but, as mentioned in the previous 
section, the analysis is more complicated due to the need to calculate (19). The derivation 
of the fast algorithm is kept concise by relying as quickly as possible on the structure of the 
underlying multiplicative group, but we need to take into consideration the cost c„ of working 
with polynomials over finite fields. In Section 4 we will give efficient methods to calculate Wq. 

The product over j in Lemma 1 can be reused and extended from the previous iteration. We 
store this product in a vector Pd — {Pd{0), . . . , Pd{b"^ — 1)) of length 6™, where 

Pdih) := {[il + l,u;^{x,,,)) = Pd-iih) (^1 + 7, c^. (^v„ , 
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for all < ^ < b™ and Po{h) = 1. Thus Pd can be calculated using the stored value for Pd-i- 
Hereby we reduce the construction cost by a factor of s at the cost of 0{¥^) memory. 

The computations of ijJa{vn(h(X) q{X) / p{X))) could be done in the initialization of the al- 
gorithm. Since u„ only depends on the negative powers of X we in fact have 

f h{X)q{X) \ ( h{X)q{X) \ f h{X)q{X) mod p{X) \ 

So, for fixed p{X), we can think of Vn as being a function from '¥i,[X]/p{X) — {w{X) e 
Fb[X] : deg(u') < n} to [0, 1). We can precompute these 6" values giving a construction cost of 
0(s6"6™Cu + 6"(cij + c„)) at a cost of 0(6™ + 6") memory. However, the cost is presumably 
dominating Cuj (most certainly so for the cJq, expressions we will derive in Corollary 1). It is 
standard practice to use a lookup table based on a generator when doing multiplications over a 
finite field. Making this change the construction cost becomes 0(s&"6'" + 6"(ct^ + c„)) at a cost 
of 0{¥^ + 26") memory (we have explicitly written the constant for clarity: there is a 0{b^) cost 
for the values of uja and a 0(6") cost for the lookup table). Here, is a lot cheaper than c^, as 
one has to multiply only by the same generator to construct the table. Wc do note however that 
the 0(6") memory cost grows exponentially with a as n = am. 

For the lookup table we made use of the fact that there exists a generator g{X) for the 
multiplicative group for which 

{¥,[X]/p{X)r := {g{Xf mod p{X) : < /3 < 6" - 1} = ¥,[X]/ip{X)) \ {0}, 

since we assumed p{X) to be irreducible over Fb[X]. For brevity we define the auxiliary function 
uj to make use of the indices w.r.t. the generator g{X): 

Lo : Zb-_i ^ [0, 1) : Lo(P -~ 6) = u{l3 - S mod (6" - 1)) 

(h{X)q{X) mod p{X) 



V p{x) 

g{X)" g{X)-' modpjX) 

P{X) 
g{X)f^-^ modp(X) 



)) 



where h{X) and q{X) are such that h{X) = g{XY modp{X) and q{X) = g{X)-^ modp(X). 

Now consider the worst-case error explicitly in terms of qd{X) as Ed{qd{X)). Then we can 
write the worst-case error iteratively in the form 

EdiqdiX)) := e6-..„((<7i(X), . . . , qd{X)),p{X)) 
h=0 

^ 6"-! 

= -1+5^ E ^^-1^ + ^^ Pd-i{h)uMh{X)qd{X)/p{X))) 

= efc™,„((gi(X), . . . , qd-i{X)),p{X)) + ^ E Pd-iW ojMHX) qd{X)/p{X))) 

h=0 

= e,^,^^{{q^{X),...,qd-i{X)),p{X)) 
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+ ^a;<,(0) + ^ J2 Pd-i{h)ujMh{X)qd{X)/p{X))), 

h=l 

where the worst-case error for the zero-dimensional rule is 0. The main computational burden 
is now hidden in calculating the last sum which we can write in terms of the auxiliary function 
u) as an extended sum: 

^ uMKX)qd{X)/p{X)))Pd-i(h)^ io{P^5)Qd-i{P) (20) 

h=l 13=0 



where S is such that qd{X) = g{X) * modp(X) and 

Qd-m ■■ 



Pd-i{g{Xf mod p{X)) if deg{g{XY mod p{X)) < m, 
otherwise. 



Let Qd = {Qd{0), Qd{V^ - 2)), then we have Qd = XlJ^^^-iPd, where 
and 

1 v{X)=g{XY uiodp{X), 



otherwise. 

Thus Qd is obtained by permuting the elements of the vector (P^, 0) £ M''""^. 

This extended sum (20), calculated for all possible choices of qd{X) = g{X)~^ modp{X) G 
{¥i,[X]/p{X))'^ , i.e., < 5 < 6" — 1, is in fact a circular convolution of length 6" — 1 

6" -2 

SdiS) := ^ c^(/3 - (5 mod (fe" - 1)) Qd-iW) 

13=0 

= ^ c^(/3)Qd_i(/3 + ,5mod(6"-l)). (21) 

p=o 

Calculating this convolution in the Fourier domain by the use of fast Fourier transforms (FFTs) 
takes time 0(&"log6"), see [18] for a general reference. We obtain a construction cost for the 
fast component-by-component algorithm using FFTs of 0(s6" log fe" + 6"(cij + £„)) using 0(6") 
memory. In other words, as n = am, the factor 5™ in the original complexity has been reduced 
to alogfe™. Asymptotically this is always faster for increasing m. 

We end this section with an overview of the complexities and their memory trade of: 



Algorithm 


Construction cost 
= sjitcration cost} -I- {initialization cost} 


Memory cost 


Straightforward 


s2&"6™(c„-l-c,) 




Cache Pd vector 


s6"6'"(c„ + c,,) 




Precalculate ui 


s6"6™-f 5"(c„-f^5„) 




Fast convolution 


s6" log6" + fe"(c^ +^^,) 


6" 



All these algorithms, except the fast convolution algorithm, have iteration times of 0{b"b"^) and 
so they will all be asymptotically slower than the fast convolution algorithm. Timings on a real 



12 



machine for 6 = 2 show the break even point to be at m = 5 for a = 2 and m — 6 for a — 3. 
The fast component-by-component algorithm based on fast convolution is given in Algorithm 2. 



Algorithm 2 Fast CBC construction of higher order polynomial lattice rules 

Input: base b a prime, number of dimensions s, number of points 6'", smoothness a > 1, and 
weights 7 = (7j)i>i 

Output: Generating vector q{X) = (qi(A), . . . ,qs{X)) & Gl^ 

Choose an irreducible polynomial p{X) G Ff,[A], deg(p) = n and n = am and generator g{X) 
Sete„^0,Q„^n-,,_.(^^^^--^J 

and u = (LU,,iv„{giXY (mod p(A))/p(X)))) 

for d = 1 to s do 

Sd = ® Qd 1 (by (fast) circular convolution) 
5 ~ argmin Sd{S) 

0<i5<6"-l 

Set qdiX) = 5(A)* (mod p(A)) 
Update/set Qd and Cd = Cd-i + — uja{0) + — Sd{6) 
end for 

return q{X) = (gi(A), . . . , q,(A)) 



4 Calculation of the worst-case error 

In this section we show how to calculate the infinite sum (18) which appears in the worst-case 
error formula from Lemma 1. In Theorem 2 we show that if x can be represented exactly with 
n digit precision in base b, then u>aix) can be computed in 0{an) operations. Following that, 
in Section 4.2, Theorem 3 will state explicit forms for general x € [0, 1). More importantly, for 
b ~ 2, Corollary 1 gives explicit forms to compute u!a{x) exactly for arbitrary x and a = 2 and 3 
using elementary computer operations. 



4.1 Technical definitions 

Before we can show how to compute LOa{x) we need to introduce some technical notation which 
will be used in the proofs in the next section. To motivate the notation we first look at uja after 
expanding the definitions of walfe(x) and ra{k), see (6) and (7), and using the non-zero digit 
expansion, (4), of A: = J2t=i ^'^S where all Kq. ^ 0: 



oo 
fc=l 

oo min(Q,#fc) #fe 



'(^) = ^faik) walfe(a;) 



= Yl n 6"^°' + '^ exp(2m Ka,^a, + l/b) Yl exp(2^i Ka,^a.+l/b) . (22) 

k—1 i—1 i—a-\-l 

fc=i:fii «a. ^ <- ' 

walfc/(x) 

Due to definition (4) we have ai > ■ ■ ■ > a^k > 0, i.e., Kai is the most significant base b 
digit of fc, etc. The second product in (22) can be seen as walfc'(a;) where k' is defined by 
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k' = k - b"-' and thus < fc' < fo'^minca,**,) . Now observe that the sum over ah 

fc > 1 can be expanded into muhiple sums over ah possible digit expansions for ah fc's with r 
digits for r > 1. That is, ^k sums for the ai together with companioning sums for the Ha- from 1 
to & — 1, i.e., 

oo oo Or-i — 1 6—1 6—1 / r \ 

Y,Ka^b-%x], (23) 




7- sums s.t. r independent sums 

oo>ai>--->ai.>0 



where G{k,x) = ra{k)walk{x). 

To simphfy notation and stress the structure in what fohows, we define the fohowing triangular 
sum operator which sums over all M > ai > ■ ■ ■ > Qr > m: 

M ai-1 a^_i-l 

T„fM(5):= E E E 5(«i,...,a,.), (24) 

ai— m+r— 1 02— m+r — 2 aj.—ra 



r sums 



and formally set the zero index sum, i.e., no sums to be taken, to be the identity mapping, 

T^'mg) :=.g. 

Define the concatenation of two such operators as putting the sums next to each other: 

M at-i-1 M' a^_i-l 

(T,^(t)T5(.-t))(.g):= E E E E 5(«i, • • ■ , a.), 

ai—7n'-\-t — l at— in' at+i— m+(r — t) — 1 ar—ni 



(r — t) sums 



i.e., having two independent ranges M > ai > - ■ • > > m' and M' > at-\-i > • • • > ar ^ ra. We 
remark that although these sums might look haggardly, the interpretation of the sum operator 
by their summation range is a natural way to reason about it as the following lemma shows. 

Lemma 2. For any AI > n > m we can split TJ^ {r) into r H~ 1 sets of two independent ranges: 

r 

Ti'{r)^Y.(^,fU{t)T:^{r-t)). 
t=o 

Proof. Applying {r) to a function 5(01, . . . , a^) can be interpreted combinatorially as having 
to distribute r objects in M — m + 1 different positions, numbered from m to M, which can 
each hold at most one object and then accumulating the result of applying the function g to this 
ensemble. It is trivial to note that we can split the range in two non-overlapping ranges and 
consider all partitions of r to distribute the objects over the two ranges. □ 

Specifically we find the following expansions of this summation operator for r = 1, 2, 3: 

To-(l)=ro"-i(l) + r-(l), (25a) 
To-(2) = ro"-'(2) + r-(l)ro"-'(l) + T;r(2), (25b) 

T-(3) = rr'(3) + r-(i)rr'(2) + r-(2)ro"-i(i) + r-(3). (25c) 

As we will apply Lemma 2 to a product function, g(ai, . . . ,ar) = <7(ai) • ■ • g{a-r) it is useful 
to obtain the following result. 
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Lemma 3. If the function g{ai, ... ,ar) is of product form gi(ai) ■■■ gr{ar) , then ^(r){g), 
with n finite, can he calculated in 0{nr). 

The proof of the lemma follows from the number of operations needed in Algorithm 3. 



Algorithm 3 Compute 5*1 = Tq ^{r){gi{ai) ■ ■ • gr{ar)) hr 0{nr) operations 
Initialize 5*1 = 0, . . . , S'r = 
for a,- = to n — r do 

Sr = Sr+ gr{ar) 

for t = 1 to r 1 do 

Sr-t — Sr-t + Sr-t+l gr-tiflr + t) 

end for 
end for 
return S'l 



At the end of this algorithm we have the post conditions: 

St = r^~\r^t+l)(^g,{a,)^ , for i = 1, . . . , r. 

With a slight modification we can calculate all values of 

St=T^-\r~t+l){^g,{a,)^, for i = 1, . . . , r. (26) 

For this we just let the outer loop run up to n — 1 and make a modification in the inner loop to 
only conditionally update the value of Sr-t as long as < n — t. The modified algorithm can 
be found in Algorithm 4. This algorithm is still 0(nr). 



Algorithm 4 Compute ah St = Tq ^(r - t + l){gt{at) ■ ■ ■ ,g,.(ar)), for t = 1, . . . , r, in 0{nr) 
Initialize S'l = 0, . . . , S'r = 
for Or = to n — 1 do 

Sf = Sy ~\~ gr(^^r^ 

for t = 1 to min(r, n — Or) — 1 do 

Sr-t = Sr-t + Sr-t+l gr-tip-r + t) 

end for 
end for 

return [Si,. . . ,Sr) 



4.2 A general algorithm for x having a fixed base h precision of n 

We now consider calculating uja{x) in base h for x G [0,1) which can be represented exactly 
with n digit precision in base b: x = (0.^i^2 • • ■ ^n)fc = X^ILi'?*^ That is, x is actually a 
rational number v/b", < v < b". This is exactly the situation that occurs in the component- 
by-component construction of Section 2.5, as the w„ function (14) exactly maps the Laurent 
series over Ff,((X^^)) to rationals t;/6" with < v < 6". 
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Theorem 2. Let a,b > 2 be integers. Then for any x = vb " with n > 1 and < v < b" , the 
value of LL>aivb~" ) can be computed in at most 0{an) operations as follows: calculate the vectors 

T{x) =(T„_i,...,Ti) 

:= (^To"-i(a-l) (^nfe-('^'+i)z(x,a.)j , . . . , To"-i(l)(&-('^°-+i)z(x, a„_i)) j , 
f{x) =(f„,...,fi) 



1=1 



which can both be computed by Algorithm 4, and set 

C=(Co,...,C„_i):-^-"*ri^) 

\ i=l / t=0,...,a-l 

C" = (Co, Ci, . . . , Cq_i) (Co, Co + Ci, . . . , Co + • • ■ + Cq,_i), 
where Co = Co = 1, i^en /or < u < 6" 

{(70:0-2 • T(i;fe-") + (C„-i - 1) + C ■ f'(w6-") i/ < w < 6", 
r=l i=l 1=1 

where a ■ b denotes the dot product and for x = (0.^i^2 ■ • ■ S,n)b we set 

z{x,a^) = \ ^il"'^^ , ^' and (3{x) = - Llog6(a;)J • 

Proof. We start from expression (22). For ease of manipulation we consider two different cases 
of the base b expansions for integer k > 0: 

1. Integers k which have between 1 and (a — 1) non-zero digits in base b: 

#k 

k = ^Ka, where 1 < #A; < a - 1. 

i=l 

2. Integers k which have a or more non-zero digits in base b: 

*k a 

A; = ^ Ka, = ^ + k' , where #fc > a and < fc' < b"" . 

i=l i=l 

As such we consider, for < a; < 1, 
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oo #fc 

^«(^)= E Y[b-^''^+^'>exp{27riKa,^a,+i/b) 

fe=l i=l 
#fc<a 



fc=l i=l 
#/c>Q 

0<fc'<f)°° 



We will now expand these outer sums as in (23), but first define 



b-l 



z{x,ai) := V exp(27ri K^^^a.+iA) =\,^ !!!^"'^^ , (27) 

kJ^I l"^ if^a^ + l^O, 

to move all the independent sums, cf. (23), into the product function. Further, denote by 
f3{x) the power of of the first non-zero digit in the base 6 expansion of x G [0, 1), then the 
sum over k' for case 2 becomes 



6"°-i 6°° if aa < /3(x) - 1, i.e., X = (0.0 0* **...) 



walj;' (x) = 



at least 



fc'=o I otherwise 

=:6°°K </3(x)-l], 

where the last hne uses Iverson notation. Introducing the sum operator (24) we obtain 

a— 1 / r \ / a \ 

w„(x) = 5]To°°(r) []6-('''+i)z(x,a,) +To°°(a) fo-^^K < /^(a;) - l\\{b-^''^+^h{x,a,) . 

r=l \i=l / \ 1=1 / 

Since our function is a product function, it is convenient to only deal with the operators, which 
then shortens the notation. 

We now deal with the two cases separately. For case l,l<r<Q; — l;We apply Lemma 2 
and manipulate the following expression 

a— 1 Q — 1 r a — 1 /a — 1 \ a — 1 

Y = Y.Y.Tn{r - t) T--\t) = E E Tn{r ~ t) T-~\t) + ^ T^{r). 

r=l r=l t=0 t=l \r=t ) r=\ 

By assumption of the n digit base b precision of x the sums do not depend on x. As we 
show next; they can be calculated off line in closed form. That means we are left to deal with 
the TJ*^^(t)((7r-t+i(ar-t+i) • ■ • gr{o.r)) for i = 1, . . . , a — 1. We can use Algorithm 3 for each of 
these terms, but as they are nested, we can use Algorithm 4 to calculate them all at once in time 
0{an) upon calculating T^~^{a - 1). The T,'^ sums are given by: 

Ct \ oo ai-l at-i — 1 

]^&-('^'+i'(&-l) =(&-!)*&-* Y E b-"'--- E 

i—1 / ai— 1 a2— n+i — 2 at—n 

oo oo oo 

= (6 - !)*&"* J2 • • • E E 

at—n 02—03 + 1 ai— a2+l 
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* 6-1 



=^-"*n^- (28) 

For case 2, ^fc > a, we can also apply Lemma 2 to obtain 



t=o 

= T--\a) + T^\l)T--\a - 1) + • • • + T^Ca - l)r«-i(l) + T^{a), 
which is applied to the function 

fe-^-ta, < P{x) - 1] ^f[^-('''+i)2(x,a,)j ■ 
The sums here become 

TrW(ri6-^'^'+'^(&-l)) =&-"*n^, fori<a, 
\i=l / i=l 

and T^{a) K < /3(a;) - 1] H - 1)^ =0. 

For X ^ the condition < fi{x) ~ 1] makes it such that T!^{a) — Q as > n (and /3(x) < n 
by assumption). The other values are the same as for case 1, and wc can use the closed 
form (28). Also here we use Algorithm 4 to calculate all the sums T^^'^{a — t) in 0{an) upon 
calculating T^~^{a). 

When a; = there is no need to consider splitting at a given n. To obtain Tf^(a) we can use 
a similar derivation as for (28) to obtain a closed form: 

(a \ oo /a — 1 \ 

oc a— 1 J 

= Er^(6-i)6-(-+^)("-)n^ 

6-1 ^ - 1 



6" - 6 J- 6» - 1 ■ 

i=\ 



Again (28) can be used to calculate the T^(r) for r = 1, . . . , a — 1. This completes the proof. □ 



4.3 Explicit forms for arbitrary x and small a 

Theorem 2 uses the fact that at most the first n digits of the coordinates of the polynomial 
lattice rule can be non-zero; it is hence not surprising that the resulting computational complexity 
depends on n. Here we take a similar approach, but explicitly look at the non-zero digits of x\ this 
will turn out to be a favorable approach in case of 6 = 2, for which we find explicit expressions 
in Corollary 1. We will use the following similar notation as was set up in the beginning of 
Section 2: Let the non-zero digits base h expansion of x = (0.^i^2 ■ • •)& ^ [Oj 1) be given by 

1=1 
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where 1 < ai < • ■ • < a^x, Ca; € {Ij • ■ • : ^ ~ !}■ In particular, we will see that the power of 6 ^ 
for the most significant digit of x, i.e., ai, plays a pivotal role. For x' = we set ai = oo and 

#a; = 0. 

Theorem 3. For x G [0, 1) with non-zero digit base b expansion 



X = 

i=l 



we have 



where 



Xl^-^.^ 1 < ai < ■•■ < Ca. e {1, 1}, 



uj2ix) = si{x) + hix), 

U}3{x) = Si{x) + S2{x) + S3{x), 



#x 
.7 = 1 



and for x ^ we have 



S2{x) := b-^ -2b-''' - (ai6-ai -&)^6-"^ 

S3(x) (6-"^+i(ai6 - d - 6 + 2) - l) ^ b"'^^ - b'^aib - ai - b + Ijb-^''' siib"' x mod 1) 

+ 6~^(ai6 - ai - 6)6-2''is2(6''^a; mod 1), 
where b'^'x mod 1 ~ h°-'x — • For x = we set S2(0) = b^^ and 53(0) = 6~^(& + 1)"^- 

Proof. We start in exactly the same way as in Theorem 2, that is, we split u)a{x) into a parts 
(cf. the a — 1 parts in case 1 plus the case 2 case): 

00 a— 1 



Ja{x) = ^ra{k) walfc(a;) = ^ 5^(2;) + Sa(a 



fe=l r=l 



where 5^(2^) contains all k with exactly r digits non-zero and Saix) contains all k with at least 
a digits non-zero. We only show the derivation of the formulae for si and S2 as examples. The 
ones for §2 and S3 can be obtained similarly. With z as in (27) we find 



i=0 
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= ib-l)J2 iib - 1) + (-1)) E ^"^'^'^ K^+i ^ 0] 

1=0 e=o 

= 1 -b^b-"^. 



Likewise for S2: 



1 

6 + 1 




E 


(^+i)z(x,£) 
















-6(6- 


-2)^6-(^'- 




00 

^0] E 




(*) 


-6(6- 


00 

- 1) E b'^'' 




00 

^0] E 






-6(6 


00 

- 1) E 




00 

= 0] E ^ 







This is a combinatorial formulation in terms of the possibilities for the digits of x. The two last 
lines, marked by (**), can be combined and interpreted as summing over all possible pairs of 
digits of X of which exactly one is non-zero. This then simplifies to two decoupled sums since a 
digit cannot be at the same time zero and non-zero: 

E^-^''-*"'^K^'+i = o]E^'''^''[^^+i^o]= E^-'^'-^'^-E^-^M E^""^- 

e'=o e=o \i'=Q j=i / \j=i 



^-E^--j (^E^--j 



The other double sum, marked by (*), can also be interpreted combinatorially: the sum is taken 
over all ordered pairs of non-zero digits of x. We can write: 

CXD 00 #X #X 

t'=0 l=t + l j=l j'=j+l 



#x \ I *^ \ #x 

E E - E ^''"^ 



Thus 



52(2;) 



6+1 

-^(^-2)^(|(|E^--j (^E^-'^^j-E^ 
20 



#x 

2a j 



i=i 



5(6-1) 



b 



□ 



The case where h equals 2 is of greatest practical importance, since in that case the matrix- 
vector product (10) over F;, to generate the nodes of the QMC rule can be calculated most 
efficiently by using the bitwise operations of the computer. Additionally 

= X when 6=2. 

By specializing the previous result to 6 = 2 we obtain the following explicit formulae. 
Corollary 1. For base b ~ 2 we obtain the following explicit results: 

u}2ix) = si{x) + S2{x), 

UJ3{x) = Si{x) + S2{x) + 53(2;), 

where 

si{x) = l-2x, S2(a:) = 1/3 - 2(1 - 

hix) = (1 - 5ti)/2 + (2 - ai)x, rs^ix) = (1 - 43i2)/18 + (5ii - l)a; - (2 - ai)x'^, 
with, for < X < 1, 

ai = - Llog2(x)J , ti:=2-''\ t2: 

and ai = 0, ti = and t2 ~ when x ~ 0. 
Proof. To obtain 53(0:) we note that, for x ^ 0, 

b^'x mod 1 = b'^'x - = x/ti ~ 1, 

and thus 



-2ai 



si{x/ti - 1) = 3 - 2x/ti, S2{x/ti - 1) = 13/3 - &x/ti + 2x^ /tl 



□ 



5 Numerical tests 

We compare the explicit construction from [5], with the CBC algorithm based on (fast) circular 
convolution presented in this paper, i.e.. Algorithm 2. From [5] we note that, to obtain higher 
order digital nets of high quality, the underlying point sets in the construction should have small 
values of t. Consequently, we use Niederreiter-Xing points generated by Pirsic's implementation, 
see [23], to obtain the digital {t\m, sd)-x\etB. In Table 1 we present a typical result for 6 = 2, 
a = 2 and s = 5 and two choices of weights 7^ = 0.9-' and 7^ = The numerical data in all 
our tests shows that the new construction produces better results. 

For reference we conclude the paper with tables showing the generating vectors and worst case 
errors of higher order polynomial lattice rules in base 2 constructed using the new algorithm. All 
polynomials are given by their canonical integer representation which is the polynomial evaluated 
at X = 6 = 2. The results can be found in Table 2 and Table 3 for a = 2 and a = 3 respectively. 
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7i = 


= 0.9^ 


ecBC 


^explicit 


l3 = 




ecBC 


^explicit 


m 


= 5 


0.9291 


1.0930 


m 


= 5 


0.028917 


0.096254 


m 


= 6 


0.4085 


0.4259 


m 


= 6 


0.009912 


0.014542 


m 


= 7 


0.1778 


0.1984 


m 


= 7 


0.003427 


0.005895 


m 


= 8 


0.0747 


0.0980 


m 


= 8 


0.001175 


0.002356 


m 


= 9 


0.0312 


0.0403 


m 


= 9 


0.000406 


0.000827 


m 


= 10 


0.0128 


0.0168 


TO 


= 10 


0.000139 


0.000290 


m 


= 11 


0.0052 


0.0071 


TO 


= 11 


0.000046 


0.000091 


TO 


= 12 


0.0020 


0.0027 


TO 


= 12 


0.000014 


0.000034 



Table 1: Comparison of the worst-case errors of CBC construction and explicit construction 





b = 


2, m = 10, a = 2: 


n = 20, p = 


1179649 






j 


1 


2 


3 


4 


5 




Qj 


453270 


920860 


324514 


394664 


106142 




e 


2.14e-6 


4.55e-5 


6.27e-4 


3.75e-3 


1.30e-2 




j 


6 


7 


8 


9 


10 






587632 


279628 


676057 


626366 


856775 




e 


3.39e-2 


7.45e-2 


1.43e-l 


2.51e-l 


4.08e-l 




b = 


2, m = 12, a = 2: n = 


24, p = 28311553 




j 




1 


2 


3 


4 


5 


1j 


2028384 


13051202 


839202 


14647583 


6874738 


e 


1.34e-7 


3.44e-6 


6.58e-5 


4.72e-4 


2.02e-3 


j 




6 


7 


8 


9 


10 




6522492 


13569662 


9821234 


10570369 


406897 


e 


6.09e-3 


1.45e-2 


2.97e-2 


5.46e-2 


9.19e-2 



Table 2: Higher order rules up to 10 dimensions for 5 = 2, 7j =0.9^ and a = 2 



6 = 2, m = 7, a = 3: n = 21, p = 2621441 



j 


1 


2 


3 


4 


5 




1492861 


1022044 


1785216 


215936 


1978368 


e 


2.02e-6 


5.24e-4 


8.20e-3 


4.05e-2 


1.22e-l 


j 


6 


7 


8 


9 


10 




1197580 


1837814 


485609 


1636853 


48810 


e 


2.82e-l 


5.54e-l 


9.80e-l 


1.60 


2.48 



b = 2, m = 8, Q = 3: n = 24, p = 28311553 



j 


1 


2 


3 


4 


5 


1j 


10844342 


2604270 


5720893 


8141702 


3831799 


e 


2.51e-7 


8.85e-5 


2.43e-3 


1.45e-2 


4.95e-2 


j 


6 


7 


8 


9 


10 




3616803 


15701694 


7750425 


2240926 


493873 


e 


1.21e-l 


2.49e-l 


4.54e-l 


7.59e-l 


1.19 



Table 3: Higher order rules up to 10 dimensions for 6 = 2, 7j = 0.9^ and a = 3 
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